import tables
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

f = tables.open_file("gamma_cone10_sample.h5", "r")

Provenance information stored as attributes¶

Versions of involved code
Corsika and simtel main parameters
Simtel filenames

f.root._v_attrs

/._v_attrs (AttributeSet), 20 attributes:
   [CLASS := 'GROUP',
    CORSIKA_ver := None,
    E_max := None,
    E_min := None,
    ImageExtractor_ver := '0.2',
    PYTABLES_FORMAT_VERSION := '2.1',
    TITLE := 'Output File',
    VERSION := '1.0',
    azimuth := 0.0,
    ctapipe_ver := '0.5.2.post312',
    particle_type := 0,
    prod_site_B_field := None,
    prod_site_alt := None,
    prod_site_array := None,
    prod_site_coord := None,
    prod_site_subarray := None,
    runlist := ['/data2/deeplearning/test/gamma_20deg_0deg_run100___cta-prod3-sct_desert-2150m-Paranal-HB9-all_cone10.simtel.gz'],
    simtel_ver := None,
    spectral_index := None,
    zenith := 1.2217304706573486]

Flat file structure: everything else is stored as tables¶

Three tables containing:
- Array information
- Telescope type information
- Event information
Plus one table per telescope containing:
- 2x 1D arrays with calibrated pixels data per image (charge, arrival time)

for node in f:
    print(node)

/ (RootGroup) 'Output File'
/Array_Info (Table(289,)) 'Table of array data'
/Event_Info (Table(100,)) 'Table of Events'
/LST (Table(57,)) 'Table of LST images'
/MSTF (Table(243,)) 'Table of MSTF images'
/MSTN (Table(207,)) 'Table of MSTN images'
/MSTS (Table(188,)) 'Table of MSTS images'
/SST1 (Table(117,)) 'Table of SST1 images'
/SSTA (Table(125,)) 'Table of SSTA images'
/SSTC (Table(123,)) 'Table of SSTC images'
/Telescope_Info (Table(7,)) 'Table of telescope data'

Array information table¶

f.root.Array_Info

/Array_Info (Table(289,)) 'Table of array data'
  description := {
  "run_array_direction": Float32Col(shape=(2,), dflt=0.0, pos=0),
  "tel_id": UInt8Col(shape=(), dflt=0, pos=1),
  "tel_type": StringCol(itemsize=8, shape=(), dflt=b'', pos=2),
  "tel_x": Float32Col(shape=(), dflt=0.0, pos=3),
  "tel_y": Float32Col(shape=(), dflt=0.0, pos=4),
  "tel_z": Float32Col(shape=(), dflt=0.0, pos=5)}
  byteorder := 'little'
  chunkshape := (2259,)

tel_types=['LST','MSTS','SSTC']

arr_table=f.root.Array_Info

fig, ax = plt.subplots(figsize=(10,10))

for tel_type in tel_types:
    tel_x = [x['tel_x'] for x in arr_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')]
    tel_y = [x['tel_y'] for x in arr_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')]
    plt.scatter(tel_x, tel_y, label=tel_type)
    
ax.legend()
ax.grid()
plt.show()

Telescope information table¶

f.root.Telescope_Info

/Telescope_Info (Table(7,)) 'Table of telescope data'
  description := {
  "num_pixels": UInt32Col(shape=(), dflt=0, pos=0),
  "optical_foclen": Float32Col(shape=(), dflt=0.0, pos=1),
  "tel_type": StringCol(itemsize=8, shape=(), dflt=b'', pos=2),
  "pixel_pos": Float32Col(shape=(2, 11328), dflt=0.0, pos=3)}
  byteorder := 'little'
  chunkshape := (2,)

Event information table¶

f.root.Event_Info

/Event_Info (Table(100,)) 'Table of Events'
  description := {
  "alt": Float32Col(shape=(), dflt=0.0, pos=0),
  "az": Float32Col(shape=(), dflt=0.0, pos=1),
  "core_x": Float32Col(shape=(), dflt=0.0, pos=2),
  "core_y": Float32Col(shape=(), dflt=0.0, pos=3),
  "event_number": UInt32Col(shape=(), dflt=0, pos=4),
  "h_first_int": Float32Col(shape=(), dflt=0.0, pos=5),
  "mc_energy": Float32Col(shape=(), dflt=0.0, pos=6),
  "particle_id": UInt8Col(shape=(), dflt=0, pos=7),
  "run_number": UInt32Col(shape=(), dflt=0, pos=8),
  "LST_indices": Int32Col(shape=(4,), dflt=0, pos=9),
  "MSTF_indices": Int32Col(shape=(25,), dflt=0, pos=10),
  "MSTN_indices": Int32Col(shape=(25,), dflt=0, pos=11),
  "MSTS_indices": Int32Col(shape=(25,), dflt=0, pos=12),
  "SST1_indices": Int32Col(shape=(70,), dflt=0, pos=13),
  "SSTA_indices": Int32Col(shape=(70,), dflt=0, pos=14),
  "SSTC_indices": Int32Col(shape=(70,), dflt=0, pos=15)}
  byteorder := 'little'
  chunkshape := (110,)

Example 1: extract some infos from a given event¶

event_index = 44
my_event = f.root.Event_Info[event_index]
print('Event number: {}'.format(my_event['event_number']))
print('Energy: {} TeV'.format(my_event['mc_energy']))
print('Alt: {} rad'.format(my_event['alt']))
print('Az: {} rad'.format(my_event['az']))
tel_types = ['LST','MSTS','SSTC']
for tel_type in tel_types:
    tel_indices = '{}_indices'.format(tel_type)
    print('{} = {}'.format(tel_indices,my_event[tel_indices]))

Event number: 296407
Energy: 0.6381028294563293 TeV
Alt: 1.2547309398651123 rad
Az: 0.0808204710483551 rad
LST_indices = [ 0 22  0  0]
MSTS_indices = [76  0 77 78 79  0 80  0 81  0  0  0  0  0 82  0 83  0  0  0  0  0  0  0  0]
SSTC_indices = [ 0 53  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]

Example 2: plot all charge and arrival time images for a given event¶

tel_types=['LST','MSTS','MSTN','MSTF','SSTC','SST1','SSTA']

tel_type_size={'LST':15,'MSTS':2,'MSTN':15,'MSTF':15,'SSTC':15,'SSTA':15,'SST1':15} #Size of the scatter plot point, for visualization purposes only

tel_table=f.root.Telescope_Info

for tel_type in tel_types:
    pos_x, pos_y = [x['pixel_pos'] for x in tel_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')][0]
    tel_ids = [x['tel_id'] for x in arr_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')]
    exec('my_images = f.root.{}'.format(tel_type))
    my_indices = my_event['{}_indices'.format(tel_type)]
    for img_index in my_indices:
        if img_index > 0:
            img_charge = my_images[img_index]['image_charge']
            img_time = my_images[img_index]['image_peak_times']
            print('Telescope type: {}, ID: {}'.format(tel_type,tel_ids[list(my_indices).index(img_index)]))
            plt.figure(figsize=(15,5))
            plt.subplot(1, 2, 1)
            plt.scatter(pos_x[:len(img_charge)], pos_y[:len(img_charge)], c=img_charge, s=tel_type_size[tel_type])
            plt.axis('square')
            plt.colorbar()
            plt.subplot(1, 2, 2)
            plt.scatter(pos_x[:len(img_time)], pos_y[:len(img_time)], c=img_time, s=tel_type_size[tel_type])
            plt.axis('square')
            plt.colorbar()
            plt.show()

Telescope type: LST, ID: 2

Telescope type: MSTS, ID: 55

Telescope type: MSTS, ID: 57

Telescope type: MSTS, ID: 58

Telescope type: MSTS, ID: 59

Telescope type: MSTS, ID: 61

Telescope type: MSTS, ID: 63

Telescope type: MSTS, ID: 69

Telescope type: MSTS, ID: 71

Telescope type: MSTN, ID: 5